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C<1 Abstract 
(N 

We present a reduction method for Wilson Dirac fermions with non-zero chemical potential 
which generates a dimensionally reduced fermion matrix. The size of the reduced fermion 
matrix is independent of the temporal lattice extent and the dependence on the chemical 
potential is factored out. As a consequence the reduced matrix allows a simple evaluation 
of the Wilson fermion determinant for any value of the chemical potential and hence the 
exact projection to the canonical partition functions. 

1 Introduction 

o 

Non-perturbative lattice calculations of Quantum Chromodynamics (QCD) at zero density 
have seen remarkable progress in recent years. However, simulations at non-zero quark or 
baryon density remain a challenge due to the occurrence of a complex phase in the fermion 
determinant at non-zero chemical potential. The fluctuation of this phase is the source of 
the notorious fermionic sign problem and obstructs the straightforward simulation of the 
theory using Monte-Carlo importance sampling. This problem limits the reliability of present 
day lattice QCD calculations at finite baryon density and makes it difficult to explore the 
^ QCD phase diagram in parameter regimes which are particularly interesting, e.g. for the 

identification of different phases of matter, the determination of phase transition lines and the 
location of possible critical endpoints. However, credible non-perturbative results at non-zero 
quark density would provide important phenomenological information, e.g. for understanding 
the structure of neutron stars or the dynamics of relativistic heavy-ion collisions. 

One approach to QCD at finite density makes use of the canonical formulation where the 
net quark (or baryon) number is held constant. This can be achieved by separating the grand- 
canonical partition function Zqc into a sum of canonical partition functions Zc(k) with a 
fixed net number k of quarks and anti-quarks. Quantities at fixed chemical potential, i.e. in 
the standard grand-canonical formulation of QCD, can then be obtained by averaging over 
the canonical partition functions. It turns out that the projection to the canonical sectors can 
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be done exactly, gauge field by gauge field, however it requires the integration of the fermion 
determinant over the whole range of imaginary chemical potential <f> = i(J,/T S [0, 2tt]. 

In the past this approach could only be made practical in connection with staggered 
fermions. For those, clever fermion matrix reduction methods were developed [1, 2, 3] that 
allow the evaluation of the determinant for any value of the chemical potential once the 
eigenvalues of the reduced fermion matrix are known. In this approach the reduction of the 
fermion matrix in size by a factor half the temporal lattice extent is the crucial ingredient 
since it reduces the complexity of the eigenvalue computation by the corresponding factor 
cubed. Unfortunately, however, for Wilson fermions so far no such reduction method was 
known despite various attempts [4, 5]. 

In this paper we present such a reduction method for Wilson fermions, i.e. we derive a 
dimensionally reduced Wilson fermion matrix whose size is independent of the temporal lattice 
extent and for which the dependence on the chemical potential is factored out. It therefore 
allows easy and exact evaluation of the determinants at any value of the chemical potential and 
hence the straightforward projection to the various canonical sectors. Applications which are 
facilitated by the reduction method for Wilson fermions presented here include the reweighting 
of ensembles to different values of the chemical potential and calculations based on canonical 
ensembles [6, 7, 8, 9]. 

The reduction of the four dimensional Wilson Dirac operator to the three dimensional 
reduced fermion matrix is very similar to the construction of the four dimensional overlap 
operator from the five dimensional domain wall fermion operator [10, 11, 12]. A similar re- 
duction method for the Wilson fermion matrix has also been proposed in [13] in the context 
of reweighting with stochastic determinants. Finally, while preparing the paper we were in- 
formed by Nakamura and Nagata [14] about their development of similar reduction techniques 
for the Wilson fermion matrix. 

The paper is organized as follows. In section 2 we briefly review the separation of the 
grand-canonical partition function of QCD into a sum of canonical partition functions with 
fixed quark or baryon number. In section 3 we present the reduction method for Wilson 
fermions which renders the computational complexity of the determinant independent of the 
temporal lattice extent and factorizes the dependence on the chemical potential. In section 
4 and 5 we discuss spectral properties of the reduced matrix and some properties of the 
projected determinants, respectively. While the results from these sections so far do not have 
a direct physical application, we would like to emphasise their potential importance for the 
development of new canonical simulation algorithms, or for the optimization of reweighting 
strategies. Finally, in section 6 we present some results from a reweighting of canonical 
ensembles, merely as a demonstration of the potential of the reduced fermion matrix approach. 

2 Canonical formulation of QCD at fixed baryon number 

The grand-canonical partition function at temperature T and chemical potential \i q for a 
single quark flavor can be defined as 




(1) 
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where M(U, p q ) denotes the Dirac operator, U collects the gauge field degrees of freedom from 
the color gauge group SU(iV c ) and S g (U) is the gauge field action. This is the commonly used 
partition function for simulating QCD thermodynamics on the lattice [15, 16, 17, 18], which 
in general, however, suffers from a strong fermionic sign problem. The same thermodynamic 
physics can also be extracted using the canonical partition function [6, 7, 19, 20, 21], although 
one should keep in mind that the physics of the two systems strictly coincide only in the 
thermodynamic limit, i.e. in the limit of infinite spatial volume. The partition function in the 
canonical approach, for a system with a net number of k quarks, can be written as 

Z c (k) = J VUe- s ^ u) det k M(U), (2) 

where the fermionic contribution is now included in the projected determinant 

det fc M(Z7) = — [ d<t> e~ ik4> det M(U, p q = i<j>T) (3) 
2?r Jo 

and one has made use of the fact that det M(U, p, q = i<f>T) enjoys a ^-periodicity in cf> [22, 23]. 
The periodicity stems from the fact that a shift in the imaginary chemical potential 4> — > <i>+jif- 
can be exactly compensated by a corresponding Z(iV c )-transformation of the underlying gauge 
field. From this periodicity it also follows that Zc{k) = for k 7^ mod N c , i.e. it vanishes 
for non-integer baryon numbers tib ^ Z, while Zc{k) = Z^(—k) follows from the evenness 
Zcc(Pq) = ZGc{—^q) due to time-reversal symmetry. 

Finally, one can relate the canonical partition functions back to the grand-canonical ones 
using the fugacity expansion 

+oo 

Zgc{^)= £ e k ^ T Z c (k), (4) 

k=— oo 

where the sum can in principle be restricted to k = mod N c , i.e. integer baryon numbers, 
following the discussion above. Furthermore, the equation also motivates the determination 
of the baryon chemical potential in the canonical approach by a definition based on the 
free energy. Essentially, the baryon chemical potential is the response of the system when 
introducing one more baryon to the system, i.e. 

M"b) = F(N C ■ (n B + 1)) - F(N C ■ n B ) , (5) 

where F(k) = — T log Zc(k) is the Helmholtz free energy of the canonical partition function. 
In a finite volume V, this definition is ambiguous due to the discreteness of the baryon number, 
however, in the thermodynamic limit it yields hb(pb) = df/dps, where p B = ub/V and 
/ = F/V are the baryon and free energy densities. Note that the baryon chemical potential 
Pb is different from the quark chemical potential p q used above. The quark chemical potential 
cannot be defined as the increase of the free energy when introducing a quark in the system 
since the free energy is infinite for systems that have fractional baryon numbers. If we need 
to compare the chemical potential used in the grand-canonical approach to the chemical 
potential measured in the canonical approach, we use p q ~ pb/N c . 
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3 Reduction technique for the Wilson fermion matrix 



In the following we consider QCD with Wilson fermions on a periodic lattice with temporal 
and spatial extent L t and L s , respectively. The quark chemical potential is (jl = a\i q where a 
is the lattice spacing. 1 The massive Wilson-Dirac operator can be written as 

M = l -T u {V u + VI) - * V* V u + m, (6) 

where V„, V£ denote the covariant forward and backward lattice derivative, T u are the Eu- 
clidean Dirac matrices and m is the bare quark mass. The chemical potential fj, couples to the 
fermion number operator and is introduced on the lattice [24] by furnishing the forward and 
backward temporal hopping terms by factors of e ±A( , respectively. More explicitly we have 



M x 



, y = (m + 4) • 5 x , y - £ [p(+k) U k (x) 5 y ^ k + P(-fc) Ut(y) 5 yx _^ 

k=l 

- {e + ^P(+4) U 4 {x) 5 y>x+i + e^P{-A)u\{y) 5^} , (7) 

where U u {x) G SU(A^ C ) are the gauge field links and P(±v) = |(l=Fr^), v = 1, . . . , 4. Further, 
for convenience we introduce P± = P(±4) for the projectors in temporal direction and will 
use this notation in the following. Choosing the spatial Euclidean Dirac matrices as 

k ? ) 

and hence hermitian, we note that the spatial part of the Wilson Dirac operator, i.e. the first 
line of eq.(7), can be written in the form 

where P++ = B is hermitian and trivial in Dirac space, 

1 3 

(B ++ ) x , y = (m + 4) ■ S x , y - - ^ {S ytX+% U k (x) + 5 x y+k ul{y)} , (10) 



k=i 



while 



C*,y = {^ x+k U k (x) - 5 x y+k Ul{y)} (11) 

k=l 

is hermitian if the Heisenberg matrices a k in eq.(8) are chosen to be hermitian. The derivations 
presented below focus on the un-improved version of the Wilson Dirac operator. However, 
the addition of the clover term can be easily accommodated: it only changes the spatial 
hopping matrix B in a way that is consistent with the structure in eq.(9) which is all that is 
needed for some of the more specific derivations in sections 3.2 and 4.1. In fact, the numerical 
experiments presented later in the paper use the clover improved version of the Dirac operator. 



1 From now on we set a =1. 
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3.1 Reduction of det M 



The Wilson fermion matrix (6) in temporal gauge with (anti-) periodic boundary conditions 
in space (time) direction for a single quark flavor with chemical potential [i can be represented 
by 

/ B P+ -P- ■ ■ e-» Lt \ 

P B l P + 



M 



P- Bo 



\ -P+U 



P+ 
B Lt-l 



(12) 



/ 



where the B^s are (4 • N c ■ L s x 4 • N c ■ L s )-matrices and represent the (spatial) Wilson 
Dirac operator on time-slice t and where we have rescaled the fermion fields such that the 
dependence on the chemical potential resides on the links connecting the last and the first 
time slice. The temporal gauge links also reside only on those links and are collected in the 
matrix U so that P+ ■ U is of the same size as Bt. For convenience we abbreviate in the 
following A~ = —IA^ ■ e~^ Lt and A + = —U ■ e +flLt and note that the matrix can be written in 
the form 



M 



( B 
1 



V 



A- 



\ /B 

P- + 



1 #L t -l / 



1 

Bi 



\A+ 



\ 



B Lt -i J 



(13) 



using the fact that P + + P_ = 1. Essentially, this splits the matrix into two parts describing 
the components of the Dirac particle moving forward and backward in time. Next we define 
the shift-projection matrix 



V 



V p 



\ 



(14) 



which leaves the forward moving part invariant and shifts the backward moving part by on 
time slice. We further note that detP = 1. Multiplying M with V from the right we find 



/ Qo(P-A- + P+ 



M -V 



Qo 

Qi Qt 
Q 2 



\ Qi t -i(p+A + + p- 



\ 



' • Qit-2 



(15) 



where we defined 



Qf = BiP T + P± 



(16) 
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and used the fact that 

P-A- + B P+ = Q (P-A-+P + ), 

P + A+ + B Lt ^ = Q+_ 1 (P + A+ + P_). 

Now we define the block diagonal matrix 

/ Qo 

Q = 



(17) 
(18) 



Qi 



\ 



Q Lt -i ) 



(19) 



and find 



/ (P_A-+P + 



M=QT X -M-V 



To 

1 Ti 
1 



V r Li _ 1 (p + A+ + p_) 



T L t -2 
1 J 



(20) 



where 



T i = (Qr)- 1 . Q +. 



(21) 



Note that the matrix M is essentially a transfer matrix describing fermions hopping forward 
and backward between time slices. We discuss this further in section 3.3. 

We can now easily calculate the determinant of the transfer matrix M using Schur com- 
plement techniques [12]. Defining 



T • Ti • . . . • T Lt 



we find 



and hence 



det [Q- 1 ■ M ■ V] = det [(P-A~ + P+) - {-l) Lt T ■ (P + A + + P_)] 



det [M] 



H detQr -det [Q^-M-V] . 



(22) 



(23) 



(24) 



i=0 



Note that the first factor det[Q] from canceling the effect of multiplication with Q~ l is inde- 
pendent of \x and cancels when we take the ratio of two determinants, e.g. with two different 
chemical potentials. Furthermore, from now on we assume Lt to be even in order to get rid 
of the inconvenient factor (— 1) *. 

In order to separate the dependence on the chemical potential /x from the gauge field de- 
pendence we first note that for arbitrary matrices A, B, C, D diagonal in Dirac space, i.e. com- 
muting with P±, we have 

(P±A + P T B)(P±C + P T D) = P±AC + P T BD . (25) 

So multiplying eq.(23) by det[P + A + P^B] with A = e~^ Lt and B = -U we obtain 

det [QT 1 -M-V] det[P+e^ Lt - PM\ = det [e~^ Lt +T-U] . (26) 
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This is the determinant of the reduced Wilson fermion matrix. Note that the gauge field 
dependence resides in T -U only and is completely separated from the dependence on [i. This 
allows now for an efficient calculation of the determinant as a function of \x for a fixed gauge 
field background. Denoting the eigenvalues of T ■ U by Aj, i = 1, . . . , 4:N C L^ we have 

4N c Ll 

det [e-» Lt +T-U]= \[ {e~» Lt + A*). (27) 

i=l 

In order to establish the equivalence between det[M] and eq.(27) we need to cancel the 
contribution from multiplying with det[P+e _/xLt — P-IA]. First we note that the matrix has 
an inverse, 

(P + e + » Lt - P_«t) • (P + e-^ Lt - PM) = 1, (28) 

hence the canceling is always possible. In fact we can calculate the determinant explicitly, 
since the matrix splits into the two orthogonal blocks oc P± (this is due to the fact that U and 
e~^ Lt are trivial in Dirac space) and the determinant is just the product of the determinants 
of the two subblocks, 

det[P+e^ L ' - PM] = det[e-^*] • det[-W] = e~' iLt ' 2NoL - (29) 

where the factor of 2 in the exponent comes from the fact that the subblock matrix spans 
over only half the Dirac indices. 



3.2 Calculation of T 

A drawback of the matrix reduction is that we need to explicitly calculate T which contains 
the inverses of of size {&N C L? S ) x (4iV c Lf). It turns out, however, that the size of the 
matrix to be inverted can be halved. To see this consider the following explicit form of QJ . 
The spatial Wilson Dirac operator Bi on time slice i inherits the structure of the full spatial 
Wilson Dirac operator B according to eq.(9). When the in are chosen to be hermitian 
it can be written as 

A Ci 



with D\ = A and C\ = d 2 . Then, in a basis where I^ is block diagonal, one has 

QT=P_ + B i P + = ( ^ J) . (31) 
We then find for the inverse of Qj 

«*">"' = ( cSr 1 !)' (32) 

so Bi needs to be inverted only in the subspace proportional to P + , i.e. only D~ l is needed. 
Further we also need to calculate det[Q^~] although these factors cancel exactly in the ratio 



2 Note that a similar argument goes through for the choice of antihermitian Dirac matrices in which case 
one has C„- = — CV 
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of determinants, e.g. for different chemical potentials. From eq.(31) we read off det[Q i ] = 
det[Dj]. Finally, for later use we also note the explicit form of Qf, 

Qt = P+ + B ,P^(l «). (33) 



3.3 Physical interpretation of the reduced matrix 

It is straightforward to give a physical interpretation of the dimensionally reduced Wilson 
Dirac fermion matrix. The equivalence of the determinants 



det M = det Q ■ det 



(34) 



establishes an equivalence (up to the bulk term det Q) between the four-dimensional Wilson 
Dirac operator M and the effective three-dimensional operator e~^ Lt ^ 2 + T -U ■ e + ^ Lt ^ 2 , similar 
to the equivalence between the five-dimensional domain wall fermion operator and the four- 
dimensional overlap Dirac operator [10, 11, 12]. The analogy becomes even more transparent 
at fi = when the reduced operator becomes 1 + T ■ U, similar to the massless overlap Dirac 
operator. The factors e ±/xLt//2 can then be understood as mass terms for the fermions and 
anti-fermions propagating forward and backward in time. 

This picture is corroborated by the matrix M in eq.(20) which is essentially a transfer 
matrix describing fermions and (anti-)fermions hopping forward and backward between time 
slices, respectively. The dynamics of these hoppings are described by the transfer matrices 
Ti which, however, are independent of /i. Fermions winding around the time direction in 
forward direction will eventually pick up a factor e +llLt (residing in A + ) for each winding 
while the anti-fermions winding around the time direction in backward direction will pick up 
corresponding factors of e~^ Lt (residing in A~). In addition, the winding fermion modes are 
then weighted by U which contains the temporal gauge field dynamics. In that sense, the 
reduced matrix is equivalent to a fermion winding number expansion. 

To complement this interpretation it is worthwile to consider alternative forms of the 
reduced matrix. In the present form the reduced matrix e~ flLt ^ 2 + T - U ■ e + ^ Lt / 2 makes the 
propagation of the fermion forward in time explicit. One can equally well emphasize the 
propagation of the anti-fermion backwards in time. This can for example be achieved by 
considering, instead of Q _1 • M ■ V, the dimensional reduction of V • M • Q _1 where Q is the 
block diagonal matrix containing Qf along the diagonal. The reduced matrix then becomes 

+ U ] . f . (35) 

where 

f = f Lt - 1 -...-f 1 'f Q with f t = Tr\ (36) 

so making the backward propagation of the anti-fermions (and their weighting with —pt instead 
of explicit. 3 

Finally, the construction of the reduced matrix presented in section 3.1 is done for gauge 
field configurations fixed to temporal gauge. However, the construction can easily be extended 



3 We note that further equivalent variants of the reduced matrix can be obtained. By considering projection 
of M with V instead of V from the left or right, one is lead to reduced matrices with modified X" or T' = J 1 ' - 1 
related to the original T by T' = T t . 
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to the generic case without gauge fixing, leading to a reduced matrix where the structure T ■ U 
becomes 

U-i 

T ■ U ■ Ti ■ Ui ■ . . . ■ T u - X ■ Z4 t -i = H Ti ■ Hi . (37) 

4=0 

Here the matrices lit now collect all the temporal gauge links at fixed time coordinate t, so 
the matrix lih t -\ is just U from before. The factors in the product of eq.(37) can be cyclically 
permuted without changing the physical content, i.e. the spectrum of the reduced matrix. 
This is due to the fact that all the cyclic permutations are related to each other by similarity 
transformations (involving the matrices Hi and Ti which have determinant one, cf. section 
4.1) while the first term in the reduced matrix is trivial in Dirac and color space. 



4 Spectral properties of the reduced matrix 



4.1 Symmetry of T • U 

In the construction of the reduced Wilson fermion matrix the crucial object is the matrix 
T -li. It turns out that this matrix has interesting properties which express themselves in 
peculiar symmetry properties of the eigenvalue spectrum. 

The first thing to note is that 

det T ■ U = 1 . (38) 

This can easily be seen from eq.(32) and eq.(33) where we read off det(Q~)~ 1 = det DJ 1 and 
detQf = det Di, respectively, and hence detTj = det [(Q^)^ 1 • Qf] = 1. 

Secondly, we note that the eigenvalues of T ■ IA come in pairs: for every eigenvalue of A, 
there is an eigenvalue A' = 1/A*. This can be seen as follows. The product Tj = (Q^)^ 1 ■ Qf 
can be LDU-decomposed with the help of eq.(32) and eq.(33): 



i 



d; 1 







Di 



1 Q 
1 



In this form it is easy to calculate its hermitian conjugated inverse, i.e. 



t 



-a 
1 



Di 




Comparing this with eq.(39) we find that 



{(Qzr 1 - Qtv 1 ]* = s - {{orr 1 - ot) ■ s- 





D7 1 



with S 



1 

-a 



and hence 



(T-U) 



-i 



t 



S-(T-U)-S- 



(39) 



(40) 



(41) 



(42) 



As a consequence the matrix T-U shares the eigenvalue spectrum with its hermitian conjugated 
inverse, that is, for each eigenvalue A G spec (T ■ li) there is another eigenvalue 1/A* € 
spec (T • li). The spectral symmetry hints at the possibility that the reduced matrix could 
be further compressed in size by a factor of two without loosing any spectral information. In 
principle this can be achieved by the projection of T ■ U to a suitable subspace, but so far 
we have not been able to construct such a projection, essentially due to the fact that T -U is 
non-normal. 



9 



Figure 1: The eigenvalue distribution for 10 arbitrary configurations drawn from canonical ensembles with 
tib — 4 (top) and ns = 11 (bottom). The red point indicates a scaled value of the spatially averaged Polyakov 
loop to show its correlation with the eigenvalue distribution. 

4.2 Eigenvalue distribution of the reduced matrix 

So far we have been concerned with purely algebraic properties of the Wilson Dirac matrix M 
and the corresponding reduced matrix T-U. In practice, what is needed are all the eigenvalues 
Aj of the latter matrix, such that the determinant of M([i) can be evaluated for any arbitrary 
value of the chemical potential according to 

4N C L 3 S 

det M(p) = det Q ■ e + ^ 2N ^ J} (e"^' + A*). (43) 

i=l 

Apart from the symmetry property A •<-)• 1/A* discussed in section 4.1 the eigenvalue spectrum 
has additional interesting features. To illustrate these we will use configurations generated 
for a previous study of non-zero baryon density systems using the canonical partition func- 
tion [25]. These 6 3 x 4 configurations are picked from Nf = 4 ensembles at a temperature close 
to the deconfining transition: T ~ 0.95T C . The parameters for the fermionic matrix will be 
set to the values used to generate the ensembles: k = 0.1371 and c sw = 1.96551 corresponding 
to a pion mass of about 700 — 800 MeV. 

In Fig. 1 we plot the eigenvalue distribution of the matrix T ■ U for a random selection of 
gauge field configurations drawn from canonical ensembles with baryon number tib = 4 (top 
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Figure 2: The eigenvalue spectrum of T ■ U for an arbitrary configuration with argP((7) ~ (middle plot), 
in the case when U = f is put by hand (left plot), or when all spatial links are set to one (right plot). 



two rows) and n# = 11 (bottom two rows), respectively. The plots represent the eigenvalue 
distribution in the complex plane and the scale is from Re, Im A £ [—600, 600]. Only the large 
magnitude eigenvalues are visible (the low magnitude members are also plotted but they are 
not visible on this scale). Note that there is a cone empty for each configuration and that 
the distributions exhibit a three lob structure which, from configuration to configuration, is 
related by Z(3)-rotations in the complex plane. An interesting observation is the fact that 
for each configuration the structure is correlated with the value of the spatially averaged 
Polyakov loop, which in the temporal gauge is given by 

p M = Jj4 trU - < 44 > 

In Fig. 1 we make this correlation explicit for each configuration by imposing the correspond- 
ing (rescaled) value of P(U) onto the eigenvalue spectrum. Of course this correlation is no 
surprise, since from the structure T-IA it is immediately clear that under a Z(3)-transformation 
of the temporal gauge fields U the eigenvalues are simply rotated by the corresponding Z(3) 
factor. On the other hand, one should also keep in mind that some of the correlations we 
observe between the determinant and the Polyakov loop might be due to the rather heavy 
quark mass [26] and the correlations might become less pronounced as we move towards lower 
quark masses. 

In order to further expose the influence of the Z(3) phase of P(U) (or rather IX) on the 
spectrum, we perform the following exercise. Instead of computing the eigenvalues of the 
original reduced matrix T -U, we calculate the spectrum of a modified reduced matrix where 
we set the temporal gauge fields U to the Z(3) phase as given by the Polyakov loop P(U). We 
show the result of this calculation for a configuration with argP({7) ~ in Fig. 2 in the left 
most plot. The only gauge field dependence is now through the spatial gauge links in T and 
we see that this dependence is responsible for the eigenvalues' variation in magnitude, while 
the phase of the eigenvalues fluctuates very weakly around zero. (Note that the spectrum 
of the reduced free Wilson Dirac operator is real and the eigenvalues A > 1 span the range 
between roughly 6 and 2000.) Of course we can also turn the argument around and put 
all the spatial links to unity, so that the gauge field dependence resides in U alone. The 
spectrum of the reduced matrix modified in this way is shown in the right most plot in Fig. 2. 
Comparing this with the original spectrum shown in the middle plot we conclude that the 
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phase variation of the eigenvalues is determined almost solely by the temporal gauge fields U 
while the spatial gauge fields only add small fluctuations to the phase. In Fig. 11 we repeat 
this exercise for configurations with &rgP(U) ~ ±2n/3. As discussed above, in this case the 
spectra are simply rotated by the corresponding Z(3) factors, and the conclusion remains the 
same. 

5 Projected determinant of the canonical partition function 

After having discussed the properties of the reduced matrix and its spectrum, the next inter- 
esting quantity to study is the determinant as a function of the (real or imaginary) chemical 
potential, and eventually also the projected determinant which subsumes the dynamics of the 
fermions in the canonical partition function. 

5.1 The determinant at non-zero chemical potential 

Having the eigenvalues A, of T ■ IA at hand it is now easy to evaluate the determinant of the 
Wilson Dirac operator for arbitrary chemical potential according to eq.(43). In Fig. 3 we 
show the logarithm of the determinant det M{n) for three configurations from a canonical 
ensemble with ub = 4. From top to bottom the configurations have argP(C7) ~ 0, +2tt/2> 
and — 27r/3. The first row shows the logarithm of |detM(/i)| normalized to detM(/x = 0), 
i.e. | log det M(/j,) | /det M(/i = 0), while the second row shows the argument argdetM(/i) 
modulo 2ir. We note significant differences between the results for configurations in the 
various triality sectors. Firstly, while in all sectors the size of the determinant ratio starts 
to grow exponentially with \fi\ beyond \fi\ > 1, configurations in the non-trivial Z(3)-sectors 
show a local minimum just below |^| ~ 1 in contrast to configurations in the trivial sector 
which show a monotonic increase with \fj,\. Secondly, the derivative of the phase w.r.t. \x stays 
roughly constant for |^| < 0.5, but depends strongly on the Z(3)-sector. 

In Fig. 4 we show the same kind of plots for three typical configurations (with argP({7) ~ 
0, +27r/3 and — 2-7r/3 from top to bottom) from a canonical ensemble with ub = 11- The 
picture qualitatively remains the same except that the features described before are accen- 
tuated. For configurations in the non-trivial Z(3)-sectors the minima of | detM(/i)| become 
slightly deeper and move to slightly larger values of \fj,\. Furthermore, for those configura- 
tions the phase of det M(^) changes more rapidly while the change of the phase in the trivial 
Z(3)-sector becomes smoother and stretches further into larger values of \fi\, possibly allowing 
reweighting to larger chemical potential. 

Note that the wild phase fluctuations observed for configurations with non-trivial Polyakov 
loop are due to the contributions from the fractional baryon number sectors. To illustrate 
this point, following the ideas presented in [23], we define a modified determinant that only 
includes the integral baryon number sectors: 



This definition is equivalent to the fugacity expansion for the determinant where we only sum 
over the sectors that have a net number of quarks divisible by 3. In figure 5 we plot the absolute 
value and phase of this modified determinant for a configuration with argP(C7) ~ 27r/3 (the 
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Figure 3: The logarithm of the determinant |detM(/x)| normalised to detM((i = 0) (left row) and the 
argument of det M{fi) (right row) for three configurations from a canonical ensemble with iib = 4. From top 
to bottom the configurations have argP((7) ~ 0, +2tv/3 and ~2n/3. 



same one used in the middle panel of figure 4). We see that the magnitude now increases 
monotonically and that the phase changes much slower with fi - this is exactly the behavior 
of configurations that have arg P(U) ~ 0. 

As we emphasized in the introduction these observations are hard to interpret physically, 
but we believe that they might play an important role for the optimization of reweighting 
strategies, or for the development of new canonical simulation algorithms. 



5.2 Calculation of the projected determinants 

Using eq.(43) we can show that the projected determinants detfc M defined in eq.(3), up to a 
multiplication with det Q, are the coefficients Cfc+fc max of the polynomial 

U(x) = 11 (x + Xi) = CkX k , (46) 
i=l fc=0 
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Figure 5: Same as the middle panel of Fig. 4 but for the modified determinant eq.(45). 



where fc max = 2N c Ng. A couple of coefficients can be computed easily: C2fc max = 1 and 
Co = Y\iK = det T • U = 1 as discussed before. All other coefficients can be calculated 
recursively. To show this, we first define the partial product: 

n n (x) = H(x + \ i ) = Y,ct ) x k - (47) 

i<n k<n 
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Figure 6: Simple model calculation: on the left the exact values and the model values are shown, 
right panel the ratio between the exact and model values are shown. 
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Clearly, IT n+ i(x) = (x + X n+ i)U n (x) and hence we have 
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(48) 



for all < k < n + 1 (we set c_l = 0). For Hq all coefficients are zero except for 



1. 



Using eq.(48) to compute c^ ri+1 - > from c^ n \ after 2k max steps we obtain the coefficients of 

n 2 fe = n. 

The recursive steps of the iteration have to be carried out using a high precision library. 
The reason for this is that the magnitude of the coefficients vary over thousands of orders of 
magnitude. We used GNU Multi-Precision library which can easily handle numbers of this 
magnitude. While the calculation takes significantly longer than when using the standard 
floating point arithmetic, the total time is small compared to the time it takes to compute 
the eigenvalues of the reduced matrix. One possible issue when using a high precision library 
is that the results look deceptively precise since the inputs are treated as high precision 
numbers when their real precision is at the level of machine precision or less. To check the 
robustness of our calculation we performed the following tests. In the first test we added 
random numbers of the order of 10~ 15 (double precision level) to the links and recomputed 
the projected determinants; the relative change in the projected determinants was at the level 
of 10~ 10 . Next, we randomly reordered the eigenvalues of the reduced matrix and repeated 
the recursive step. We find that the relative change in the coefficients is of the order of 10~ 9 . 
We conclude that our procedure is robust and it produces accurate results. 

Before we conclude this section, we want to point out that the A «-)■ 1/A* symmetry of 
the reduced matrix, together with the fact that det T ■ U = 1 can be used to show that 
c fc m ax+fc = c fc -fc- Since det Q is real, this insures that that det^M = (det_fcM)* a fact 
easily derived from the definition of the projected determinant and the reality of detM(^). 



5.3 Size distribution of the projected determinant 

One of the interesting aspects of this calculation is the fact that the magnitude of the projected 
determinants varies over many orders of magnitude as we change the quark number sector. 
In this section we will show that the bulk of this variance can be captured by a simple 
combinatorial argument. However, this only captures part of the variance: in order to describe 
the variance accurately we need to take into account the complex phase of the eigenvalues of 
the reduced matrix. 
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Figure 7: Comparison between true polynomial coefficients and the ones where the eigenvalue phase is 
quenched. 

In the main, the bulk of the magnitude is given by the average magnitude of the eigen- 
values and combinatorics. To illustrate this point, in Fig. 6 we compare the value of the 
coefficients of the polynomial U(x) for a given configuration with the coefficients of the poly- 
nomial n'(x) = (x + A) fcmax (x + 1/A) fcmax , where A is the geometric mean of the magnitude of 

the large eigenvalues A = ( 11[A|>1 IM ) 

For any given configurations the eigenvalues vary both in phase and magnitude; in the 
comparison above we quenched both fluctuations. While this approximation captures a good 
part of the variation of magnitude, there is still quite a discrepancy left. To trace the source 
of discrepancy we compare the coefficients of II(x) with a polynomial where each eigenvalue 
is replaced with its magnitude: II" (x) = n«( x + I'M)- The results of this comparison are 
presented in Fig. 7: we see that the discrepancy is similar to the one above. This shows that 
the source of discrepancy is the complex phase fluctuation that is disregarded here. This 
conclusion is also supported by the fact that the coefficients of II" (x) are larger in magnitude 
than the coefficients of n(x); this is due to cancelations produced by phase fluctuations. 




Figure 8: Comparison between true polynomial coefficients and the ones where the eigenvalue phase is drawn 
from a random Z(3) distribution. Note the difference in scale in the right plot as compared to Fig. 7 and 6. 

To prove that phase fluctuations are responsible for the discrepancy we model the rough 
features of the phase distribution. The phase distribution is illustrated in Fig. 1. One obvious 
feature is that the eigenvalues are concentrated around the Z(3) axes. This makes sense: if 
the eigenvalues were exactly along the Z(3) axes, the projected determinants that have zero 
triality {k = 3ng) will be real whereas the ones with non-zero triality will have a fluctuating 
phase and they will vanish when averaged over gauge configurations, as expected. 
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To check that phase fluctuations are responsible for reducing the magnitude of the co- 
efficients we compare the coefficients of U(x) with U'"(x) = Yl^x + Zj|Aj|), where Z{ is a a 
random Z(3) phase with one constraint: we enforce the (A,l/A*) pairing. The results are 
presented in figure 8: we see that the results agree much better now. It is clear that the phase 
fluctuations play an important role in determining the coefficients. Phase fluctuations reduce 
the value of the coefficients significantly: from figure 7 we see that for intermediate k values 
they reduce the magnitude by about 200 orders of magnitude. From figure 8 we see that this 
is mostly due to the Z(3) nature of the fluctuations (we also tested random 7L{n) phases with 
n 6 {2, 4, 5, 6} but they fail to produce the same agreement). 



6 Phase fluctuations in canonical ensembles 



The reduced Wilson fermion matrix constructed and discussed in the previous sections allows 
for several intriguing applications. Presumably among the most interesting are the direct 
simulation of canonical ensembles and the reweighting to different fermion numbers, or values 
of the chemical potential. In the following we briefly discuss these two applications and 
present some results on the phase fluctuations encountered in direct simulations of canonical 
ensembles, merely to illustrate the capabilities of the reduced fermion matrix approach. 

In order to simulate the canonical partition function eq.(2) by Monte-Carlo techniques 
one needs the integrand to be real and positive. Since in general this is not guaranteed, the 
approach so far [7] has been to ensure positivity by fiat, i.e. to generate an ensemble using 
the weight W\k\(U) tx |RedetfcM(£/")|, while the phase 

= det k M(U) 

a[U) ~ w w (u) i4yj 

is introduced when computing the observable 

where denotes the average with regard to the generated ensemble based on the 
measure. In practice, to evaluate the partition function numerically, the continuous Fourier 
transform in eq.(3) has so far either been replaced by a discrete Fourier transform [7] or 
by a more sophisticated approximation [27, 25] and the so introduced bias needs a careful 
treatment. With the reduced Wilson fermion matrix these approximations have become 
obsolete, simply because the Fourier transform can now be evaluated exactly. 

Still, the quenching of the phase of the integration measure, det&M(?7) — > W\u, introduces 
a systematic error which one needs to control. A measure of how severe the fluctuations of 
the phase are, is provided by the expectation value of a(U) in eq.(49) in a given ensemble. It 
turns out that the sign fluctuations for the canonical ensembles seem to be under good control 
and one might wonder how generic this feature is. One way to estimate the reliability of the 
simulations is to reweigh results generated at one value of k to other values k' and to check 
consistency between the reweighted results and the ones obtained from direct simulations. For 
the reweighting from one canonical ensemble Zc(k) to another Zc(k') the relevant quantity 
is 

det fc ,M(EQ 

O|*|-*(E0= w {U) ■ (51) 
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For this definition of the reweighting factor we see from the second column of Table 1 that 
its magnitude changes very fast as we move away from the original ensemble, in this case 
the one with rig = 4. However, the value of the factor changes for all configurations in a 
similar manner and the average is still comfortably away from zero in terms of its variance 
(cf. third and fourth column) even though its magnitude is dramatically changed. (Note 
that this dramatic change is related to the variation of detfc M with k over many orders of 
magnitude.) What we mean by that is that its magnitude, as compared with its variance, 
is larger than two or more. In principle it is when this ratio becomes close to one that the 
reweighting in the equation above will run into numerical difficulties. We see from Table 1 
that for baryon numbers as large as n' B = 16 the average is still twice the variance. 
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Table 1: The real part of the average reweighting factor eq.(51). The imaginary part vanishes by symmetry. 
This factor is computed based on the ensemble generated with ub = 4 (marked with a star above). 



In Fig. 9 we plot the ratio between the variance and the average of the real phase factor 
reweighted from the canonical ensembles with ng = 4 and 11, as a function of the reweighted 
baryon number. Note that in both cases the average reweighting factor is well behaved over 
a large range of reweighted baryon numbers. 
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Figure 9: The real part of the average phase factor: the ratio between its variance and its mean value. 



Encouraged by this result one might also try to reweight the chemical potential as a 
function of the baryon number following the definition in eq.(5). The chemical potential, as 
defined there, can be shown to be 

f s/rr i z c(3(n B + 1)) , / det 3( n B +i)M \ 

HB(n B )/T = -log — — — = -log( — — — ) (52) 

Zc{3n B ) \ det 3 n fl M / 3ns 

<det 3(rifl+1) M/^|3 nB |> 

= ~ lQ g Ta I tuTFvu \ • ( 53 ) 

(det 3nB M/W l3nBl ) m 
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In Fig. 10 we show the reweighted chemical potential as defined by eq.(5), for the same 
ensembles as before, ng = 4 (left plot) and ng = 11 (right plot). Despite the absence of a 
sign problem in the reweighting factor, it is evident that the reweighting fails. The ensemble 
riB = 4 in the confined phase misses to describe even the neigbouring ensemble with only 
slightly different baryon number ub = 5. The ensemble ub = 11, on the other hand, is in 
the deconfined phase and seems to be able to describe other deconfined phases with different 
baryon numbers. Not surprisingly, however, it completely fails to describe ensembles with 
tib < 8, and hence the phase transition, simply because it does not contain any information 
about the confined phase. 





Figure 10: The baryon chemical potential reweighted from the ensemble with ub = 4 (left plot) and ub = 11 
(right plot). Empty circles are the results of the direct calculation. 



The next interesting question would now be to reweight an ensemble in the mixed phase, 
e.g. with ns = 7. However, the reweighted results turn out to be too noisy, indicating that the 
information gathered in the canonical ensembles is simply not enough to allow for a reliable 
reweighting. One has to keep in mind that the canonical ensembles used here only contain 
about 1500 configurations each. The next obvious step is then to employ a multi-ensemble 
reweighting, combining the information from all the ensembles into the reweighting. What 
we find is that for n' B = 4 the most important contributions indeed come from the ensembles 
with small ns and that, as we move towards larger n' B , higher ensembles come into play as 
expected. However, even the multi-ensemble reweighting does not seem capable to reproduce 
the S-shape behaviour typical for a first order phase transition, and we conclude that the 
reweighting suffers from a severe overlap problem. Nevertheless, the example impressively 
demonstrates the ease with which such calculations can be achieved using the reduced Wilson 
fermion matrix. 



7 Conclusions 



We have presented a reduction method for Wilson Dirac fermions which generates a dimen- 
sionally reduced fermion matrix. The size of the reduced matrix is independent of the temporal 
lattice extent. Moreover, the dependence of the matrix on the chemical potential factors out 
and reduces to a simple multiplicative factor. This allows to evaluate the Wilson fermion de- 
terminant for any value of the chemical potential, once the eigenvalues of the reduced matrix 
are calculated, and hence allows to perform the exact projection of the determinant to the 
canonical sectors with fixed fermion number. 
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The reduced fermion matrix presented here facilitates various interesting applications, 
for example the reweighting of ensembles to arbitrary values of the chemical potential or 
the fermion number. So far, this has only been possible for staggered fermions. Another 
application is the direct simulation of canonical ensembles and this is now possible without 
any bias from inexact projections to the canonical sectors. Since the size of the reduced 
matrix is independent of the temporal lattice extent, such calculations can in principle be 
done at arbitrarily low temperatures, barring possible sign problems. 

The reduced fermion matrix has some interesting properties like the spectral symmetry 
A o 1/A*, a simple behavior of the spectrum under Z(A r c )-transformations and the correlation 
of the spectrum with the Polyakov loop. We believe that such properties may be important 
for the development of more efficient canonical simulation algorithms, or for the optimisation 
of reweighting strategies. 

As a first test we applied the reduction method to a set of canonical ensembles and 
determined the phase fluctuations of the Wilson fermion determinant at non-zero chemical 
potential, or non-zero fermion number, using standard reweighting techniques. It turns out 
that for the ensembles considered here, the overlap problem inherent in all reweighting meth- 
ods introduces a systematic bias that forbids reliable calculations, e.g. using multi-ensemble 
reweighting. On the other hand, the phase fluctuations seem to be rather well controlled and 
it will be interesting to see whether this is a generic feature of canonical partition functions. 
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